<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Spectral factorization using Kolmogorov 1939 approach.</title>
<link rel="canonical" href="/Users/mcgrant/Repos/CVX/examples/filter_design/html/spectral_fact.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Spectral factorization using Kolmogorov 1939 approach.</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
Text output
&nbsp;&nbsp;&nbsp;&nbsp;
Plots
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% (code follows pp. 232-233, Signal Analysis, by A. Papoulis)</span>
<span class="comment">%</span>
<span class="comment">% Computes the minimum-phase impulse response which satisfies</span>
<span class="comment">% given auto-correlation.</span>
<span class="comment">%</span>
<span class="comment">% Input:</span>
<span class="comment">%   r: top-half of the auto-correlation coefficients</span>
<span class="comment">%      starts from 0th element to end of the auto-corelation</span>
<span class="comment">%      should be passed in as a column vector</span>
<span class="comment">% Output</span>
<span class="comment">%   h: impulse response that gives the desired auto-correlation</span>

<span class="keyword">function</span> h = spectral_fact(r)

<span class="comment">% length of the impulse response sequence</span>
n = length(r);

<span class="comment">% over-sampling factor</span>
mult_factor = 100;        <span class="comment">% should have mult_factor*(n) &gt;&gt; n</span>
m = mult_factor*n;

<span class="comment">% computation method:</span>
<span class="comment">% H(exp(jTw)) = alpha(w) + j*phi(w)</span>
<span class="comment">% where alpha(w) = 1/2*ln(R(w)) and phi(w) = Hilbert_trans(alpha(w))</span>

<span class="comment">% compute 1/2*ln(R(w))</span>
w = 2*pi*[0:m-1]/m;
R = [ ones(m,1) 2*cos(kron(w',[1:n-1])) ]*r;
alpha = 1/2*log(R);

<span class="comment">% find the Hilbert transform</span>
alphatmp = fft(alpha);
alphatmp(floor(m/2)+1:m) = -alphatmp(floor(m/2)+1:m);
alphatmp(1) = 0;
alphatmp(floor(m/2)+1) = 0;
phi = real(ifft(j*alphatmp));

<span class="comment">% now retrieve the original sampling</span>
index  = find(rem([0:m-1],mult_factor)==0);
alpha1 = alpha(index);
phi1   = phi(index);

<span class="comment">% compute the impulse response (inverse Fourier transform)</span>
h = real(ifft(exp(alpha1+j*phi1),n));
</pre>
</div>
</body>
</html>